The No-U-Turn Sampler 



The No-U-Turn Sampler: Adaptively Setting Path Lengths 

in Hamiltonian Monte Carlo 

Matthew D. Hoffman mdhoffma@cs.princeton.edu 



Department of Statistics 

Columbia University 

New York, NY 10027, USA 

Andrew Gelman 

Departments of Statistics and Political Science 

Columbia University 

New York, NY 10027, USA 



GELMAN@STAT.COLUMBIA.EDU 



Abstract 

Hamiltonian Monte Carlo (HMC) is a Markov chain Monte Carlo (MCMC) algorithm that 
avoids the random walk behavior and sensitivity to correlated parameters that plague many 
MCMC methods by taking a series of steps informed by first-order gradient information. 
These features allow it to converge to high-dimensional target distributions much more 
quickly than simpler methods such as random walk Metropolis or Gibbs sampling. However, 
HMC's performance is highly sensitive to two user-specified parameters: a step size e 
and a desired number of steps L. In particular, if L is too small then the algorithm 
exhibits undesirable random walk behavior, while if L is too large the algorithm wastes 
computation. We introduce the No-U-Turn Sampler (NUTS), an extension to HMC that 
eliminates the need to set a number of steps L. NUTS uses a recursive algorithm to build 
a set of likely candidate points that spans a wide swath of the target distribution, stopping 
automatically when it starts to double back and retrace its steps. Empirically, NUTS 
perform at least as efficiently as and sometimes more efficiently than a well tuned standard 
HMC method, without requiring user intervention or costly tuning runs. We also derive a 
method for adapting the step size parameter e on the fly based on primal-dual averaging. 
NUTS can thus be used with no hand-tuning at all. NUTS is also suitable for applications 
such as BUGS-style automatic inference engines that require efficient "turnkey" sampling 
algorithms. 

Keywords: Markov chain Monte Carlo, Hamiltonian Monte Carlo, Bayesian inference, 
adaptive Monte Carlo, dual averaging. 



1. Introduction 

Hierarchical Bayesian models are a mainstay of the machine learning and statistics com- 
munities. Exact posterior inference in such models is rarely tractable, however, and so 
researchers and practitioners must usually resort to approximate statistical inference meth- 



ods. Deterministic approximate inference algorithms (for example, those reviewed by Wain 



wrig ht and Jorda^ ( |2008D ) can be efficient, but introduce bias and can be difficult to apply 



to some models. Rather than computing a deterministic approximation to a target poste- 
rior (or other) distribution, Markov chain Monte Carlo (MCMC) methods offer schemes for 
drawing a series of correlated samples that will converge in distribution to the target distri- 
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bution (Neal, 1993). MCMC methods are sometimes less efficient than their deterministic 



counterparts, but are more generally applicable and are asymptotically unbiased. 

Not all MCMC algorithms are created equal. For complicated models with many param- 



eters, simple methods such as random- walk Metropolis (Metropolis et al. , 1953) and Gibbs 



sampling (Geman and Geman, 1984) may require an unacceptably long time to converge 



to the target distribution. This is in large part due to the tendency of these methods to 



explore parameter space via inefficient random walks (Neal, 1993). When model parameters 



are continuous rather than discrete, Hamiltonian Monte Carlo (HMC), also known as hybrid 
Monte Carlo, is able to suppress such random walk behavior by means of a clever auxiliary 
variable scheme that transforms the problem of sampling from a target distribution into the 



problem of simulating Hamiltonian dynamics (Neal, 2011). The cost of HMC per indepen- 



dent sample from a target distribution of dimension D is roughly 0{D^^^), which stands in 
sharp contrast with the 0{D^) cost of random-walk Metropolis (Creutz 1988). 

HMC's increased efficiency comes at a price. First, HMC requires the gradient of the 
log-posterior. Computing the gradient for a complex model is at best tedious and at worst 
impossible, but this requirement can be made less onerous by using automatic differentiation 
(Griewank and Walther, 2008). Second, HMC requires that the user specify at least two 



parameters: a step size e and a number of steps L for which to run a simulated Hamiltonian 
system. A poor choice of either of these parameters will result in a dramatic drop in HMC's 
efficiency. Methods from the adaptive MCMC literature (see Andrieu and Thoms (2008|) for 
a review) can be used to tune e on the fly, but setting L typically requires one or more costly 
tuning runs, as well as the expertise to interpret the results of those tuning runs. This hurdle 
limits the more widespread use of HMC, and makes it challenging to incorporate HMC into 



a general-purpose inference engine such as BUGS (Gilks and Spiegelhalter 1992), JAGS 



( htt p : / / mcmc-j ags . sour c eforge . net ) , Infer.NET (Minka et al. ), HBC (Daume HI, 2007), or 



PyMC (Patil et al., 2010). 



The main contribution of this paper is the No- U- Turn Sampler (NUTS), an MCMC 
algorithm that closely resembles HMC, but eliminates the need to choose the problematic 
number-of-steps parameter L. We also provide a new dual averaging (Nesterov 2009) 



scheme for automatically tuning the step size parameter e in both HMC and NUTS, making 
it possible to run NUTS with no hand-tuning at all. We will show that the tuning-free 
version of NUTS samples as efficiently as (and sometimes more efficiently than) HMC, even 
ignoring the cost of finding optimal tuning parameters for HMC. Thus, NUTS brings the 
efficiency of HMC to users (and generic inference systems) that are unable or disinclined to 
spend time tweaking an MCMC algorithm. 



2. Hamiltonian Monte Carlo 



In Hamiltonian Monte Carlo (HMC) ( |Neal[|201lHl993t|Duane et al.[|1987[ ), we introduce an 
auxiliary momentum variable for each model variable 9d- In the usual implementation, 
these momentum variables are drawn independently from the standard normal distribution, 
yielding the (unnormalized) joint density 



p{6, r) oc exp{£(^) 



r}. 



(1) 



2 



The No-U-Turn Sampler 



Algorithm 1 Hamiltonian Monte Carlo 



Given 0°, e, L, £,M: 
for m = 1 to M do 

Sample ~ AA(0,/). 

Set 6"^ ^ ^ ^ f ^ r". 

for i = 1 to L do 

Set 6,r ^ Leapfrog(0, f, e). 

end for 

With probability a = min il, ,^p|7(eii^-!^_'I^?^o} ^ set 0" 



end for 

function Leapfrog(^, r, e) 
Set f ^r + (e/2)V0£(0). 
Set 6i ^ 6* + ef. 
Set f + (e/2)Ve£(^). 
return 0,f. 



where C is the logarithm of the joint density of the variables of interest (up to a normalizing 
constant) and x ■ y denotes the inner product of the vectors x and y. We can interpret this 
augmented model in physical terms as a fictitious Hamiltonian system where 9 denotes a 
particle's position in D-dimensional space, denotes the momentum of that particle in 
the dth dimension, £ is a position-dependent negative potential energy function, ■ r is 
the kinetic energy of the particle, and \ogp{6, r) is the negative energy of the particle. We 
can simulate the evolution over time of the Hamiltonian dynamics of this system via the 
"leapfrog" integrator, which proceeds according to the updates 



r 



t+e/2 ^^t^ (e/2)Ve£(0*); 9'+' = 9' + er'+'/^; r'+' = r'+'/^ + {e/2)VeC{9'+'), (2) 



where r* and 0* denote the values of the momentum and position variables r and 9 at time 
t and Vg denotes the gradient with respect to 9. Since the update for each coordinate 
depends only on the other coordinates, the leapfrog updates are volume-preserving — that 
is, the volume of a region remains unchanged after mapping each point in that region to a 
new point via the leapfrog integrator. 

A standard procedure for drawing M samples via Hamiltonian Monte Carlo is described 
in Algorithm [TJ / denotes the identity matrix and A/'(/x, S) denotes a multivariate normal 
distribution with mean /i and covariance matrix S. For each sample m, we first resample 
the momentum variables from a standard multivariate normal, which can be inetpreted as 
a Gibbs sampling update. We then apply L leapfrog updates to the position and momen- 
tum variables 9 and r, generating a proposal position-momentum pair 0,f. We propose 
setting 9"^ = 9 and = — f, and accept or reject this proposal according to the Metropo- 



lis algorithm (Metropolis et al. 1953). This is a valid Metropolis proposal because it is 



time-reversible and the leapfrog integrator is volume-preserving; using an algorithm for 
simulating Hamiltonian dynamics that did not preserve volume would seriously complicate 
the computation of the Metropolis acceptance probability. The negation of f in the pro- 
posal is theoretically necessary to produce time-reversibility, but can be omitted in practice 
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if one is only interested in sampling from p{6). The algorithm's original name, "Hybrid 
Monte Carlo," refers to the hybrid approach of alternating between updating 9 and r via 
Hamiltonian simulation and updating r via Gibbs sampling. 

The term log , on which the acceptance probability a depends, is the negative 
change in energy of the simulated Hamiltonian system from time to time eL. If we could 
simulate the Hamiltonian dynamics exactly, then a would always be 1, since energy is con- 
served in Hamiltonian systems. The error introduced by using a discrete-time simulation 
depends on the step size parameter e — specifically, the change in energy | log ^j§7y | is pro- 
portional to for large L, or if L = 1 (Leimkuhler and Reich, 2004). In theory the 
error can grow without bound as a function of L, but in practice it typically does not when 
using the leapfrog discretization. This allows us to run HMC with many leapfrog steps, 
generating proposals for 9 that have high probability of acceptance even though they are 
distant from the previous sample. 

The performance of HMC depends strongly on choosing suitable values for e and L. If 
e is too large, then the simulation will be inaccurate and yield low acceptance rates. If e 
is too small, then computation will be wasted taking many small steps. If L is too small, 
then successive samples will be close to one another, resulting in undesirable random walk 
behavior and slow mixing. If L is too large, then HMC will generate trajectories that loop 
back and retrace their steps. This is doubly wasteful, since work is being done to bring the 
proposal 9 closer to the initial position 9^~^ . Worse, if L is chosen so that the parameters 
jump from one side of the space to the other each iteration, then the Markov chain may 
not even be ergodic (Neal, 2011). More realistically, an unfortunate choice of L may result 
in a chain that is ergodic but slow to move between regions of low and high density. 



3. Eliminating the Need to Hand-Tune HMC 



HMC is a powerful algorithm, but its usefulness is limited by the need to tune the step size 
parameter e and number of steps L. Tuning these parameters for any particular problem re- 
quires some expertise, and usually one or more preliminary runs. Selecting L is particularly 
problematic; it is difficult to find a simple metric for when a trajectory is too short, too long, 
or "just right," and so practitioners commonly rely on heuristics based on autocorrelation 



statistics from preliminary runs (Neal, 2011). 



Below, we present the No- U- Turn Sampler (NUTS), an extension of HMC that eliminates 
the need to specify a fixed value of L. In section [3^2] we present schemes for setting e based 



on the dual averaging algorithm of Nesterov ( 2009 ) 



3.1 No-U-Turn Hamiltonian Monte Carlo 

Our first goal is to devise an MCMC sampler that retains HMC's ability to suppress random 
walk behavior without the need to set the number L of leapfrog steps that the algorithm 
takes to generate a proposal. We need some criterion to tell us when we have simulated 
the dynamics for "long enough," i.e., when running the simulation for more steps would 
no longer increase the distance between the proposal 9 and the initial value of 9. We use 
a convenient criterion based on the dot product between f (the current momentum) and 
9 — 9 (the vector from our initial position to our current position), which is the derivative 
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Figure 1: Example of building a binary tree via repeated doubling. Each doubling proceeds 
by choosing a direction (forwards or backwards in time) uniformly at random, 
then simulating Hamiltonian dynamics for 2^ leapfrog steps in that direction, 
where j is the number of previous doublings (and the height of the binary tree). 
The figures at top show a trajectory in two dimensions (with corresponding binary 
tree in dashed lines) as it evolves over four doublings, and the figures below show 
the evolution of the binary tree. In this example, the directions chosen were 
forward (light orange node), backward (yellow nodes), backward (blue nodes), 
and forward (green nodes). 



with respect to time (in the Hamiltonian system) of half the squared distance between the 
initial position 9 and the current position 9: 

In other words, if we were to run the simulation for an infinitesimal amount of additional 
time, then this quantity is proportional to the progress we would make away from our 
starting point 9. 

This suggests an algorithm in which one runs leapfrog steps until the quantity in equation 
[3] becomes less than 0; such an approach would simulate the system's dynamics until the 
proposal location 9 started to move back towards 9. Unfortunately this algorithm does 
not guarantee time reversibility, and is therefore not guaranteed to converge to the correct 
distribution. NUTS overcomes this issue by means of a recursive algorithm reminiscent of 



the doubling procedure devised by Neal (2003) for slice sampling 



NUTS begins by introducing a slice variable u with conditional distribution p(u\9,r) = 
Uniform(?/; [0, exp{/^(0) — • r}]), which renders the conditional distribution p{9,r\u) = 
Uniform(0, r; {9', r'\ exp{£(^) — ■ r} > u}). This slice sampling step is not strictly neces- 
sary, but it simplifies both the derivation and the implementation of NUTS. In addition to 
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Figure 2: Example of a trajectory generated during one iteration of NUTS. The blue ellipse 
is a contour of the target distribution, the black open circles are the positions 9 
traced out by the leapfrog integrator and associated with elements of the set of 
visited states B, the black solid circle is the starting position, the red solid circles 
are positions associated with states that must be excluded from the set C of 
possible next samples because their joint probability is below the slice variable n, 
and the positions with a red "x" through them correspond to states that must be 
excluded from C to satisfy detailed balance. The blue arrow is the vector from the 
positions associated with the leftmost to the rightmost leaf nodes in the rightmost 
height-3 subtree, and the magenta arrow is the (normalized) momentum vector 
at the final state in the trajectory. The doubling process stops here, since the 
blue and magenta arrows make an angle of more than 90 degrees. The crossed- 
out nodes with a red "x" are in the right half-tree, and must be ignored when 
choosing the next sample. 



being more complicated, the analogous algorithm that eliminates the slice variable seems 
empirically to be slightly less efficient than the algorithm presented in this paper. 
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At a high level, after resampling u\9, r, NUTS uses the leapfrog integrator to trace out a 
path forwards and backwards in fictitious time, first running forwards or backwards 1 step, 
then forwards or backwards 2 steps, then forwards or backwards 4 steps, etc. This doubling 
process implicitly builds a balanced binary tree whose leaf nodes correspond to position- 
momentum states, as illustrated in Figure[T} The doubling is halted when the subtrajectory 
from the leftmost to the rightmost nodes of any balanced subtree of the overall binary tree 
starts to double back on itself (i.e., the fictional particle starts to make a "U-turn"). At 
this point NUTS stops the simulation and samples from among the set of points computed 
during the simulation, taking care to preserve detailed balance. Figure [2] illustrates an 
example of a trajectory computed during an iteration of NUTS. 

Pseudocode implementing a efficient version of NUTS is provided in Algorithm |3j A 
detailed derivation follows below, along with a simplified version of the algorithm that 
motivates and builds intuition about Algorithm [s] (but uses much more memory and makes 
smaller jumps). 



3.1.1 Derivation of simplified NUTS algorithm 



NUTS further augments the model p{6, r) oc exp{£(0) — ^r-r} with a slice variable u (Neal 



2003). The joint probability of 0, r, and u is 

p{e, r, u) oc l[u £ [0, exp{£(6') - • r}]], (4) 

where ![•] is 1 if the expression in brackets is true and if it is false. The (unnormalized) 
marginal probability of and r (integrating over u) is 

p{6, r) oc exp{£(6') — ■ r}, (5) 

as in standard HMC. The conditional probabilities p{u\9, r) and p{6, r\u) are each uniform, 
so long as the condition u < exp{£(6') — |r • r} is satisfied. 

We also add a finite set C of candidate position-momentum states and another finite set 
B ^ C to the model. B will be the set of all position-momentum states that the leapfrog 
integrator traces out during a given NUTS iteration, and C will be the subset of those 
states to which we can transition without violating detailed balance. B will be built up by 
randomly taking forward and backward leapfrog steps, and C will selected deterministically 
from B. The random procedure for building B and C given 9, r, u, and e will define a 
conditional distribution p{B,C\6, r, u, e), upon which we place the following conditions: 

C.l: All elements of C must be chosen in a way that preserves volume. That is, any 
deterministic transformations of 0, r used to add a state 9\ r' to C must have a Jacobian 
with unit determinant. 

C.2: p{{9,r) e C|0,r,u,e) = 1. 

C.3: p{u < exp{£{9') - \r' ■ r'}|(0',r') G C) = 1. 

C.4: If (0,r) G C and {9' ,r') G C then for any B, p{B,C\d,r,u, e) = p{B,C\9' ,r' ,u, e). 

C.l ensures that p{9, r\{9, r) G C) oc p{9, r), i.e. if we restrict our attention to the elements 
of C then we can treat the unnormalized probability density of a particular element of C as 
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an unnormalized probability mass. C.2 says that the current state 9,r must be included in 
C. C.3 requires that any state in C be in the slice defined by u, i.e., that any state {9' ,r') £ C 
must have equal (and positive) conditional probability density p{9',r'\u). C.4 states that B 
and C must have equal probability of being selected regardless of the current state 9, r as 
long as {9,r) £ C (which it must be by C.2). 

Deferring for the moment the question of how to construct and sample from a distribu- 
tion p{B,C\9, r, u, e) that satisfies these conditions, we will now show that the the following 
procedure leaves the joint distribution p{9,r,u,B,C\e) invariant: 

1. sample r ~ AA(0, /), 

2. sample u ~ Uniform([0, exp{£(^*) — ■ r}]), 

3. sample B, C from their conditional distribution p{B, C\9^, r, u, e), 

4. sample 6l*+\r ~ T(6i*,r,C), 

where T{9' ,r'\9,r,C) is a transition kernel that leaves the uniform distribution over C in- 
variant, i.e., T must satisfy 

^ E n«'./i«..c).ffi'^ (0) 

' ' {e,r)ec ' ' 

for any 9',r'. The notation 0*+^,r ~ T(9^,r,C) denotes that we are resampling r in a way 
that depends on its current value. 

Steps 1, 2, and 3 resample r, u, B, and C from their conditional joint distribution given 
9^, and therefore together constitute a valid Gibbs sampling update. Step 4 is valid because 
the joint distribution of 9 and r given u, B, C, and e is uniform on the elements of C: 

p{9, r\u, B, C, e) oc p{B,C\9, r, n, e)p{9, r\u) 

oc p{B,C\9,r,u,e)I[u < exp{C{9) - ■ r}] (7) 
ocI[{9,r) G C]. 

Condition C.l allows us to treat the unnormalized conditional density p{9,r\u) oc I[u < 
exjp{C{9) — • r}] as an unnormalized conditional probability mass function. Conditions 
C.2 and C.4 ensure that p{B,C\9,r,u,e) oc I[{9,r) £ C] because by C.2 {9,r) must be in C, 
and by C.4 for any B, C pair p(B,C\9, r, u, e) is constant as a function of 9 and r as long as 
{9, r) £ C. Condition C.3 ensures that {9,r) £ C ^ u < exp{C{9) — • r} (so the p{9, r\u,e) 
term is redundant). Thus, equation [7] implies that the joint distribution of 9 and r given u 
and C is uniform on the elements of C, and we are free to choose a new 0*+^, r*"^^ from any 
transition kernel that leaves this uniform distribution on C invariant. 

We now turn our attention to the specific form for p{B,C\9,r,u, e) used by NUTS. 
Conceptually, the generative process for building B proceeds by repeatedly doubling the 
size of a binary tree whose leaves correspond to position-momentum states. These states 
will constitute the elements of B. The initial tree has a single node corresponding to the 
initial state. Doubling proceeds by choosing a random direction Vj ~ Uniform({ — 1, 1}) and 
taking 2^ leapfrog steps of size vje (i.e., forwards in fictional time if Vj = 1 and backwards in 
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fictional time if Vj = — 1), where j is the current height of the tree. (The initial single-node 
tree is defined to have height 0.) For example, if vj = 1, the left half of the new tree is the 
old tree and the right half of the new tree is a balanced binary tree of height j whose leaf 
nodes correspond to the 2^ position-momentum states visited by the new leapfrog trajectory. 
This doubling process is illustrated in Figure [T] Given the initial state 9, r and the step size 
e, there are 2^ possible trees of height j that can be built according to this procedure, each 
of which is equally likely. Conversely, the probability of reconstructing a particular tree of 
height j starting from any leaf node of that tree is 2~^ regardless of which leaf node we 
start from. 

We cannot keep expanding the tree forever, of course. We want to continue expanding B 
until one end of the trajectory we are simulating makes a "U-turn" and begins to loop back 
towards another position on the trajectory. At that point continuing the simulation is likely 
to be wasteful, since the trajectory will retrace its steps and visit locations in parameter 
space close to those we have already visited. We also want to stop expanding B if the 
error in the simulation becomes extremely large, indicating that any states discovered by 
continuing the simulation longer are likely to have astronomically low probability. (This 
may happen if we use a step size e that is too large, or if the target distribution includes 
hard constraints that make the log-density C go to — oo in some regions.) 

The second rule is easy to formalize — we simply stop doubling if the tree includes a leaf 
node whose state 6, r satisfies 

/:(^) -^r-r -log n<-Amax (8) 

for some nonnegative Amax- We recommend setting Amax to a large value like 1000 so 
that it does not interfere with the algorithm so long as the simulation is even moderately 
accurate. 

We must be careful when defining the first rule so that we can build a sampler that 
neither violates detailed balance nor introduces excessive computational overhead. To de- 
termine whether to stop doubling the tree at height j, NUTS considers the 2^ — 1 balanced 
binary subtrees of the height-j tree that have height greater than 0. NUTS stops the dou- 
bling process when for one of these subtrees the states 9~,r~ and 9~^,r^ associated with 
the leftmost and rightmost leaves of that subtree satisfies 

{9+ -9~)-r' <0 or {9+ - 9') ■ r+ < 0. (9) 

That is, we stop if continuing the simulation an infinitesimal amount either forward or back- 
ward in time would reduce the distance between the position vectors 9~ and 9~^. Evaluating 
the condition in equation [o] for each balanced subtree of a tree of height j requires 2-'+^ — 2 
inner products, which is comparable to the number of inner products required by the 2^ — 1 
leapfrog steps needed to compute the trajectory. Except for very simple models with very 
little data, the cost of these inner products is usually negligible compared to the cost of 
computing gradients. 

This doubling process defines a distribution p{B\9, r, u, e). We now define a deterministic 
process for deciding which elements of B go in the candidate set C, taking care to satisfy 
conditions C.1~C.4 on p(B,C\9,r,u,€) laid out above. C.l is automatically satisfied, since 
leapfrog steps are volume preserving and any element of C must be within some number 
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of leapfrog steps of every other element of C C.2 is satisfied as long as we include the 
initial state 9,r in C, and C.3 is satisfied if we exclude any element 6\r' of B for which 
exp{£(0') — ^r' ■ r'} < u. To satisfy condition C.4, we must ensure that p{B,C\9,r,u,e) = 
p{B,C\9' ,r' ,u,e) for any {9',r') G C. For any start state {9',r') E B, there is at most one 
series of directions {vq, . . . ,Vj} for which the doubling process will reproduce B, so as long 
as we choose C deterministically given B either p{B,C\0' ,r' ,u,e) = 2"-' = p{B,C\9,r,u,€) 
or p{B,C\9' ,r' ,u, e) = 0. Thus, condition C.4 will be satisfied as long as we exclude from 
C any state 9',r' that could not have generated B. The only way such a state can arise is 
if starting from 9',r' results in the stopping conditions in equations [s] or [o] being satisfied 
before the entire tree has been built, causing the doubling process to stop too early. There 
are two cases to consider: 

1. The doubling procedure was stopped because either equation [8] or equation [9] was 
satisfied by a state or subtree added during the final doubling iteration. In this case 
we must exclude from C any element of B that was added during this final doubling 
iteration, since starting the doubling process from one of these would lead to a stopping 
condition being satisfied before the full tree corresponding to B has been built. 

2. The doubling procedure was stopped because equation [9] was satisfied for the leftmost 
and rightmost leaves of the full tree corresponding to B. In this case no stopping 
condition was met by any state or subtree until B had been completed, and condition 
C.4 is automatically satisfied. 

Algorithm [2] shows how to construct C incrementally while building B. After resam- 
pling the initial momentum and slice variables, it uses a recursive procedure resembling a 
depth-first search that eliminates the need to explicitly store the tree used by the doubling 
procedure. The BuildTree() function takes as input an initial position 9 and momentum r, 
a slice variable u, a direction v G { — 1, 1}, a depth j, and a step size e. It takes 2^ leapfrog 
steps of size ve (i.e. forwards in time if v = 1 and backwards in time if t; = —1), and returns 

1. the backwardmost and forwardmost position-momentum states 9~,r~ and 9^,r~^ 
among the 2^ new states visited; 

2. a set C of position-momentum states containing each newly visited state 9',r' for 
which exp{£(6'') — ^r' • r'} > u; and 

3. an indicator variable s; s = indicates that a stopping criterion was met by some state 
or subtree of the subtree corresponding to the 2^ new states visited by BuildTree(). 

At the top level, NUTS repeatedly calls BuildTree() to double the number of points that 
have been considered until either BuildTree() returns s = (in which case doubling stops 
and the new set C' that was just returned must be ignored) or equation [9] is satisfied for 
the new backwardmost and forwardmost position- momentum states 9~,r~^ and 9~^,r^ yet 
considered (in which case doubling stops but we can use the new set C). Finally, we select 
the next position and momentum 9"^, r uniformly at random from C, the union of all of the 
valid sets C that have been returned, which clearly leaves the uniform distribution over C 
invariant. 
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Algorithm 2 Naive No-U-Turn Sampler 

Given 9°, e, £, M: 
for m = 1 to M do 

Resample r° ^ 7V(0, /). 

Resample u ^ Uniforni([0, exp{C{6"^~^ — ^r" • r"}]) 

Initialize 0- = 9""-^, 9+ = 9"^-^, r~ = r^, r+ = r^, j = 0, C = {(6l™-i, rO)}, s = 1. 
while s — \ do 

Choose a direction Vj ^ Uniforni({— 1, 1}). 

if Vj = —\ then 

9~ , — , — ,C', s' -i— BuildTree(6'^, , u, Vj, j, e). 

else 

-,-,9+,r+,C',s' ^ BuildTree(6l+,r+,w,t;j,j,e). 
end if 

if s' = 1 then 
C^CUC. 
end if 

s ^ s%{9+ -9-)-r- > 0]l[{9+ - ff") • r+ > 0]. 
end while 

Sample 9"^,r uniformly at random from C. 
end for 

function BuildTree(0, r, u, v, j, e) 
if j = then 

Base case — take one leapfrog step in the direction v. 
9' ,r' -tr- Leapfrog(6', r, ve). 
rr^j W^r')} ifu<eM^O')-lr' -r'} 
^ \ else 
s' ^ l[u < exp{A,nax + ^9') - \r' ■ r'}]. 
return 9\r\9' ,r' ,C' ,s' . 
else 

Recursion — build the left and right subtrees. 
9-,r-,9+,r+,C',s' ^ BuildTree(6l, r, u, w, j - l,e). 
if V — —1 then 

9-,r-,-,~,C",s" ^ BuildTree(6'-,r-,u,w, j - l,e). 
else 

-,-,9+,r+X'\s" ^ BuildTree(6'+, r+,u,u,j - l,e). 
end if 

s' ^ s's"I[(6i+ - 9-) ■ r- > 0]I[{9+ - 9-) • r+ > 0]. 
C'-^C'UC". 

return 9^ ,r^ , 9^, r^,C' , s' . 
end if 



To summarize, Algorithm [2] defines a transition kernel that leaves p{9, r, u, B,C\e) invari- 
ant, and therefore leaves the target distribution p{6) oc exp{£(^)} invariant. It does so by 
resampling the momentum and slice variables r and u, simulating a Hamiltonian trajectory 
forwards and backwards in time until that trajectory either begins retracing its steps or 
encounters a state with very low probability, carefully selecting a subset C of the states 
encountered on that trajectory that lie within the slice defined by the slice variable u, and 
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finally choosing the next position and momentum variables 9"^ and r uniformly at random 
from C. Figure [2] shows an example of a trajectory generated by an iteration of NUTS where 
equation |9] is satisfied by the height-3 subtree at the end of the trajectory. Below, we will 
introduce some improvements to algorithm [2] that boost the algorithm's memory efficiency 
and allow it to make larger jumps on average. 

3.1.2 Efficient NUTS 

Algorithm [2] requires 2-' — 1 evaluations of C{6) and its gradient (where j is the number 
of times BuildTree() is called), and 0{2^) additional operations to determine when to stop 
doubling. In practice, for all but the smallest problems the cost of computing £ and its 
gradient still dominates the overhead costs, so the computational cost of algorithm [2] per 
leapfrog step is comparable to that of a standard HMC algorithm. However, Algorithm 
[2] also requires that we store 2-' position and momentum vectors, which may require an 
unacceptably large amount of memory. Furthermore, there are alternative transition kernels 
that satisfy detailed balance with respect to the uniform distribution on C that produce 
larger jumps on average than simple uniform sampling. Finally, if a stopping criterion 
is satisfied in the middle of the final doubling iteration then there is no point in wasting 
computation to build up a set C that will never be used. 

The third issue is easily addressed — if we break out of the recursion as soon as we 
encounter a zero value for the stop indicator s then the correctness of the algorithm is 
unaffected and we save some computation. We can address the second issue by using a more 
sophisticated transition kernel to move from one state {6, r) G C to another state {9' , r') G C 
while leaving the uniform distribution over C invariant. This kernel admits a memory- 
efficient implementation that only requires that we store 0{j) position and momentum 
vectors, rather than 0{2^). 

Consider the transition kernel 

T(z.KC) = <; + = if|c-i<|c-| ' 



where w and w' are shorthands for position-momentum states [0, r), C"°"' and C"'"^ are disjoint 
subsets of C such that (^n°»'ij(joid _ ^^^^ ^ ^ ^oid_ English, T proposes a move from C"'"^ 
to a random state in C"""™ and accepts the move with probability ^^om^ ■ This is equivalent 

to a Metropolis-Hastings kernel with proposal distribution q{w' ,0°^"^' \w,C°^'^,C'"''") oc 
I[w' € C'"'^]I[C°'''' = C"<="]]I[C"<=-' = C°'''], and it is straightforward to show that it satisfies 
detailed balance with respect to the uniform distribution on C, i.e. 

p{w\C)T{w'\w,C) =p{w'\C)T{w\w',C), (11) 

and that T therefore leaves the uniform distribution over C invariant. If we let C'^"^ be 
the (possibly empty) set of elements added to C during the final iteration of the doubling 
(i.e. those returned by the final call to BuildTree() and C"'"* be the older elements of C, 
then we can replace the uniform sampling of C at the end of Algorithm [2] with a draw 
from T{6^ ,r* ,C) and leave the uniform distribution on C invariant. In fact, we can apply T 
after every doubling, proposing a move to each new half-tree in turn. Doing so leaves the 
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uniform distribution on each partially built C invariant, and therefore does no harm to the 
invariance of the uniform distribution on the fully built set C. Repeatedly applying T in 
this way increases the probability that we will jump to a state 9^~^^ far from the initial state 
9*; considering the process in reverse, it is as though we first tried to jump to the other 
side of C, then if that failed tried to make a more modest jump, and so on. This transition 



kernel is thus akin to delayed-rejection MCMC methods (Tierney and Mira, 1999), but in 
this setting we can avoid the usual costs associated with evaluating new proposals. 

The transition kernel above still requires that we be able to sample uniformly from the 
set C returned by BuildTree(), which may contain as many as 2^~^ elements. In fact, we 
can sample from C without maintaining the full set C in memory by exploiting the binary 
tree structure in Figure [TJ Consider a subtree of the tree explored in a call to BuildTree(), 
and let Cgubtree denote the set of its leaf states that are in C: we can factorize the probability 
that a state {9,r) £ Csubtree will be chosen uniformly at random from C' as 

P(g,rn = rj^= '^7;r'ir ^ I (12) 

I'-' I I'-' I I '-'subtree I 

= p{{9,r) £ CsnhtTee\C)p{9,r\{9,r) £ Csubtrec, C). 

That is, p{9, r\C') is the product of the probability of choosing some node from the subtree 
multiplied by the probability of choosing 9,r uniformly at random from Cgubtree- We use 
this observation to sample from C' incrementally as we build up the tree. Each subtree 
above the bottom layer is built of two smaller subtrees. For each of these smaller subtrees, 
we sample a 9,r pair from p{9, r\ (9, r) G Csubtree) to represent that subtree. We then choose 
between these two pairs, giving the pair representing each subtree weight proportional to 
how many elements of C are in that subtree. This continues until we have completed the 
subtree associated with C and we have returned a sample 9' from C and an integer weight 
n' encoding the size of C, which is all we need to apply T. This procedure only requires that 
we store 0{j) position and momentum vectors in memory, rather than 0(2-'), and requires 
that we generate 0{2^) extra random numbers (a cost that again is usually very small 
compared with the 2-^ — 1 gradient computations needed to run the leapfrog algorithm). 

Algorithm [3] implements all of the above improvements in pseudocode. Matlab code im- 
plementing the algorithm is also available at http : //www . cs . princeton . edu/~iiidhof f ma[ 
and a C++ implementation will also be available as part of the soon-to-be-released Stan 
inference package. 

3.2 Adaptively Tuning e 

Having addressed the issue of how to choose the number of steps L, we now turn our 
attention to the step size parameter e. To set e for both NUTS and HMC, we propose using 



stochastic optimization with vanishing adaptation ( Andrieu and Thoms 2008 ) , specifically 



an adaptation of the primal-dual algorithm of Nesterov (2009) 



Perhaps the most commonly used vanishing adaptation algorithm in MCMC is the 



stochastic approximation method of Robbins and Monro ( 1951 ). Suppose we have a statistic 



Ht that describes some aspect of the behavior of an MCMC algorithm at iteration t > 1, 
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Algorithm 3 Efficient No-U-Turn Sampler 

Given 9°, e, C, M: 
for TO = 1 to M do 
Resample r° ~ AA(0,/). 

Rcsample u ^ Uniform([0, exp{£(0'"~^ — ■ r°}]) 

Initialize O" = 9+ = 6''"-\ r" = r°,r+ =r°,j = 0,6'™ = 6'™-\n= l,s = 1. 

while s = 1 do 

Choose a direction Vj ^ Uniform({— 1, 1}). 

if Vj = — 1 then 

6~ ,r~ ,—,—,6' ,n' ,s' BuildTree{9~ ,r~ ,u,Vj,j, e). 

else 

-,-,9+,r+,9',n',s' ^ BuildTree(6'+, r+, w, -y^, j, e). 
end if 

if s' = 1 then 

With probabihty min{l, set 6*"' ^ 9'. 
end if 
n ^ n + n'. 

s ^ s'I[(6i+ - 6*-) • r- > 0]I[(6'+ - 61") • r+ > 0]. 

J <^ i + 1- 

end while 
end for 

function BuildTree(^, r, u, v, j, e) 
if j = then 

Base case — take one leapfrog step in the direction v. 
9' ,r' ^ Leapfrog(6', r, ve). 
n' ^ l[u < cx^p{C{9') - \r' ■ r'}]. 
s' ^ \{u < exp{A,„ax + Ci9') - \r' ■ r'}]. 
return 9' ,r' ,9' ,r' ,9' ,n' , s' . 
else 

Recursion implicitly build the left and right subtrees. 
9-,r-,9+,r+,9',n',s' ■t- BuildTree(6», r,u,v,j-l,e). 
if s' = 1 then 

if V = —1 then 

0~.r-,-,-,9",n",s" <(- BuildTree(6l-,r-,u,w, j - l,e). 

else 

-,-,9+,r+,9",n",s" <r- BuildTree(6l+, r+, u, v,i - 1, e). 
end if 

With probability ^^7^, set 9' ^ 9". 
s' ^ s"I[{9+ - 9-) ■ r- > 0]I[(6»+ - 9-) • r+ > 0] 
n' ^n' + n" 
end if 

return 9~ ,r~ ,9~^ ,r~^ ,9' ,n' , s' . 
end if 



and define its expectation h{x) as 

1 ^ 

h{x) = Et[Ht\x] = lim - J^mtlx], (13) 
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where x G M is a tunable parameter to the MCMC algorithm. For example, if at is the 
Metropolis acceptance probability for iteration t, we might define Ht = 6 — at, where 6 is 
the desired average acceptance probability. If /i is a nondecreasing function of x and a few 



other conditions such as boundedness of the iterates xt are met (see Andrieu and Thoms 



(|2008j) for details), the update 

xt+i ^ Xt- r]tHt (14) 

is guaranteed to cause h{xt) to converge to as long as the step size schedule defined by rjt 
satisfies the conditions 

J2vt = oo; ^T/j2<cx). (15) 
t t 

These conditions are satisfied by schedules of the form rjt = t~'^ for k G (0.5, 1]. As long 
as the per-iteration impact of the adaptation goes to (as it will if r]t = t~'^ and k > 0) 
the asymptotic behavior of the sampler is unchanged. That said, in practice x often gets 
"close enough" to an optimal value well before the step size r/ has gotten close enough to 
to avoid disturbing the Markov chain's stationary distribution. A common practice is 
therefore to adapt any tunable MCMC parameters during the burn-in phase, and freeze the 



tunable parameters afterwards (e.g., (Gelman et al. , 2004)). 



Dual averaging: The optimal values of the parameters to an MCMC algorithm dur- 
ing the burn-in phase and the stationary phase are often quite different. Ideally those 
parameters would therefore adapt quickly as we shift from the sampler's initial, transient 
regime to its stationary regime. However, the diminishing step sizes of Robbins-Monro give 
disproportionate weight to the early iterations, which is the opposite of what we want. 



Similar issues motivate the dual averaging scheme of Nesterov (2009), an algorithm 
for nonsmooth and stochastic convex optimization. Since solving an unconstrained convex 
optimization problem is equivalent to finding a zero of a nondecreasing function (i.e., the 
(sub)gradient of the cost function), it is straightforward to adapt dual averaging to the 
problem of MCMC adaptation by replacing stochastic gradients with the statistics Ht- 
Again assuming that we want to find a setting of a parameter x G M such that h[x) = 
¥.t[Ht\x\ = 0, we can apply the updates 

xt+i ^ n Tvrr -^5 ^ vtxt+i + (1 - 'nt)xt, (16) 

where /i is a freely chosen point that the iterates xt are shrunk towards, 7 > is a free 
parameter that controls the amount of shrinkage towards ;U, to ^ is a free parameter that 
stabilizes the initial iterations of the algorithm, rjt = t^'^ is a step size schedule obeying the 
conditions in equation [TS], and we define xi = xi. As in Robbins-Monro, the per-iteration 
impact of these updates on x goes to as t goes to infinity. Specifically, for large t we have 

xt+i-xt = 0{-Htt-^-''), (17) 

which clearly goes to as long as the statistic Ht is bounded. The sequence of averaged 
iterates xt is guaranteed to converge to a value such that h{xt) converges to 0. 

The update scheme in equation [16] is slightly more elaborate than the update scheme 



of Nesterov (2009), which implicitly has to = and k = 1. Introducing these parameters 
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addresses issues that are more important in MCMC adaptation than in more conventional 
stochastic convex optimization settings. Setting to > improves the stability of the algo- 
rithm in early iterations, which prevents us from wasting computation by trying out extreme 
values. This is particularly important for NUTS, and for HMC when simulation lengths are 
specified in terms of the overall simulation length eL instead of a fixed number of steps L. 
In both of these cases, lower values of e result in more work being done per sample, so we 
want to avoid casually trying out extremely low values of e. Setting the parameter k < 1 
allows us to give higher weight to more recent iterates and more quickly forget the iterates 
produced during the early burn-in stages. The benefits of introducing these parameters are 
less apparent in the settings originally considered by Nesterov, where the cost of a stochastic 
gradient computation is assumed to be constant and the stochastic gradients are assumed 
to be drawn i.i.d. given the parameter x. 

Allowing to > and n G (0.5, 1] does not affect the asymptotic convergence of the dual 
averaging algorithm. For any n G (0.5,1], xt will eventually converge to the same value 
i Ei=i ^t- We can rewrite the term ^^q^ as T^^toj 0{Vi), which is the 

only feature needed to guarantee convergence. 

We used the values 7 = 0.05, to = 10, and k = 0.75 for all our experiments. We arrived 
at these values by trying a few settings for each parameter by hand with NUTS and HMC 
(with simulation lengths specified in terms of eL) on the stochastic volatility model described 
below and choosing a value for each parameter that seemed to produce reasonable behavior. 
Better results might be obtained with further tweaking, but these default parameters seem 
to work consistently weU for both NUTS and HMC for all of the models that we tested. It 
is entirely possible that these parameter settings may not work as well for other sampling 
algorithms or for H statistics other than the ones described below. 

Setting e in HMC: In HMC we want to find a value for the step size e that is neither 
too small (which would waste computation by taking needlessly tiny steps) nor too large 
(which would waste computation by causing high rejection rates). A standard approach is 
to tune e so that HMC's average Metropolis acceptance probability is equal to some value 
6. Indeed, it has been shown that (under fairly strong assumptions) the optimal value of e 
for a given simulation length eL is the one that produces an average Metropolis acceptance 



probability of approximately 0.65 (Beskos et al. 2010 Neal, 2011). For HMC, we define a 
criterion /i^'^'-'(e) so that 

where 0* and f* are the proposed position and momentum at the tth iteration of the Markov 
chain, 9^~^ and r*'*^ are the initial position and (resampled) momentum for the tth iteration 
of the Markov chain, if^MC jg acceptance probability of this tth HMC proposal and 
^HMC jg ^j^g expected average acceptance probability of the chain in equilibrium for a fixed 
e. Assuming that h^^^ is nonincreasing as a function of e, we can apply the updates in 



equation 16 with Ht = 5 — ffHMC ^^^^ ^ = loge to coerce h^^'~" = 5 for any 6 G (0, 1). 



Setting e in NUTS: Since there is no single accept/reject step in NUTS we must define 
an alternative statistic to Metropolis acceptance probabihty. For each iteration we define 
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Algorithm 4 Heuristic for choosing an initial value of e 
f unct ion FindReasonableEpsilon( ) 
hiitializG e = 1, r ~ 7V(0, /). 
Set d',r' ^ Leapfrog(6', r, e). 

«^2l[^>0.5]-l. 

-hile {%^y>^- do 
e ^ 2"e. 

Set 0', r' -s- Leapfrog(0, r, e). 
end while 
return e. 



Algorithm 5 Hamiltonian Monte Carlo with Dual Averaging 
Given 6", 6, A, £, M, M^'^'^p*: 

Set £0 = FindReasonableEpsilon(6'), /I log(lOeo), eo 1, ffo = 0, 7 = 0.05, to = 10, k = 0.75. 
for m = 1 to M do 
Reample r° ~7V(0,/). 

Set 0™ ^ e"'-^,e^ e"'-\r ^ r°,L„, = max{l,Round(A/e„_i)}. 
for i = 1 to Lm do 

Set 9,r i~ Leapfrog(0, f, e„i-i)- 
end for 

With probability a = min {l, ..pffgf '^i^jf^.o} } , set 0™ ^ 0, ^ -f . 
if m < A/adapt tjjgjj 

Set i?™ = (1 - + T^i^ - «)■ 

Set log e„i = /i - ^H„i, log e,„ m"'' log + (1 - m"'") log e,„_i. 
else 

Set 67^ = 6j\,/adapt . 

end if 
end for 



the statistic i/^UTS j^-g expectation when the chain has reached equilibrium as 

where Bf^^^ is the set of all states explored during the final doubling of iteration t of the 
Markov chain and 9*~^ and r^'^ are the initial position and (resampled) momentum for the 
tth iteration of the Markov chain. ff'^UTS ^^^^ understood as the average acceptance 
probability that HMC would give to the position-momentum states explored during the final 
doubling iteration. As above, assuming that ij'^UTS jg nonincreasing in e, we can apply the 



updates in equation 16 with Ht ^ 6 — i^^UTS g^^^^ = loge to coerce /j^UTS _ ^ ^^j. g^^y 

5e (0,1). 

Finding a good initial value of e: The dual averaging scheme outlined above should 
work for any initial value ei and any setting of the shrinkage target ^. However, convergence 
will be faster if we start from a reasonable setting of these parameters. We recommend 
choosing an initial value ei according to the simple heuristic described in Algorithm |4j In 
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Algorithm 6 No-U-Turn Sampler with Dual Averaging 

Given 0°, 5, r.M.M"'*"^*: 

Set eo = FindReasonableEpsilon(^), // = log(lOeo), eo = 1, -ffo = 0, 7 = 0.05, to = 10, k = 0.75. 
for m — 1 to M do 
Sample r° ~ A/'(0,/). 

Resamplc u ~ Umform([0, cxp{/;(6''"-^ - ^r" ■ r°}]) 

Initialize = e"'-\ 9+ =9""-^, r" = r°, r+ = r°, j = 0, 6''" = n = 1, s = 1. 

while s = 1 do 

Choose a direction ~ Uniform({— 1, 1}). 

if Vj = —1 then 

9~,r~,-, -, 9', n', s', a, na BuildTree(0~ , r~ , m, Vj,j, em-i9"'~^ ,r''). 
else 

-, -, 9+,r+,9', n', s' , a, ^ BuildTree(6i+, r+, w, Vj,j, e^-i, 6»"'"\ r°). 
end if 

if s' = 1 then 

With probabihty min{l, ^}, set 6>"' 6>'. 
end if 
n •«— n + n'. 

s ^ s'I[i9+ - 9-) ■ r- > 0]1[{9+ - 9-) ■ r+ > 0]. 

end while 

if m < M'^^^'P' then 

Set H^={l- ^) + ^{6 - ^). 

Set loge™ = ^ - ^-ffmjlogem = ?n~'^logem + (1 - m"") logem-i. 
else 

Set 
end if 
end for 

function BuildTree(6i, r, u, v, j, e, 9°, r°) 
if j = then 

Base case — take one leapfrog step in the direction v. 

9' , r' -ir- Leapfrog(S, r, ve). 

n' ^ I[u < exp{£(6l') - fr' • r'}]. 

s' ^ I[u < exp{A^ax + C{e') - \r' ■ r'}]. 

return 9\ r' , 9', r' , 9', n', s', min{l, exp{£(6'') - \r' ■ r' - C{9°) + \r° ■ r°}}, 1. 
else 

Recursion — implicitly build the left and right subtrees. 
9-,r-,9+,r+,9',n',s',a',n'„ ^ BuildTree(6i, r, u, w, j - l,e,9°,r°). 
if s' = 1 then 
if w = —1 then 

9-,r-,-,-,9",n",s",a",n'^ <f- BuildTree(6i-,r-,M,v,i - l,e,6»°,r°). 
else 

-, -, 9+,r+,9",n", s" , a", < ^ BuildTree(6l+, r+, it, v, j - 1, e, g", r"). 
end if 

With probability set 9' ^ 9". 

Set q' <— a' + a" , n'a <— + n^. 
s' ^ s"I[{9+ - 9-) ■ r- > 0]I[(6»+ - 9-) ■r+ >0] 
n' <-n' + n" 
end if 

return 9~ ,r~ ,9'^ ,r'^ , 9' , n' , s', a' , n'^. 
end if 
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English, this heuristic repeatedly doubles or halves the value of ei until the acceptance 
probability of the Langevin proposal with step size ei crosses 0.5. The resulting value of ei 
will typically be small enough to produce reasonably accurate simulations but large enough 
to avoid wasting large amounts of computation. We recommend setting ^ = log(lOei), since 
this gives the dual averaging algorithm a preference for testing values of e that are larger 
than the initial value ei. Large values of e cost less to evaluate than small values of e, and 
so erring on the side of trying large values can save computation. 

Algorithms [5] and [g] show how to implement HMC (with simulation length specified in 
terms of eL rather than L) and NUTS while incorporating the dual averaging algorithm 
derived in this section, with the above initialization scheme. Algorithm [5] requires as input 
a target simulation length A ~ eL, a target mean acceptance probability (5, and a num- 
ber of iterations M^'^^^^ after which to stop the adaptation. Algorithm [g] requires only a 
target mean acceptance probability 5 and a number of iterations M''^^^^'^ . Matlab code im- 



plementing both algorithms can be found at http://www.cs.princeton.eciu/~mdhoffma 



and C++ implementations will be available as part of the Stan inference package. 



4. Empirical Evaluation 

In this section we examine the effectiveness of the dual averaging algorithm outlined in 



section 3.2, examine what values of the target 5 in the dual averaging algorithm yield 
efficient samplers, and compare the efficiency of NUTS and HMC. 

For each target distribution, we ran HMC (as implemented in algorithmjs]) and NUTS (as 
implemented in algorithm [6]) with four target distributions for 2000 iterations, allowing the 



step size e to adapt via the dual averaging updates described in section [3^ for the first 1000 
iterations. In all experiments the dual averaging parameters were set to 7 = 0.05, to = 10, 
and At = 0.75. We evaluated HMC with 10 logarithmically spaced target simulation lengths 
A per target distribution. For each target distribution the largest value of A that we tested 
was 40 times the smallest value of A that we tested, meaning that each successive A is 
40^/^ ~ 1.5 times larger than the previous A. We tried 15 evenly spaced values of the dual 
averaging target 5 between 0.25 and 0.95 for NUTS and 8 evenly spaced values of the dual 
averaging target 5 between 0.25 and 0.95 for HMC. For each sampler-simulation length-(5- 
target distribution combination we ran 10 iterations with different random seeds. In total, 
we ran 3,200 experiments with HMC and 600 experiments with NUTS. 

We measure the efficiency of each algorithm in terms of effective sample size (ESS) 
normalized by the number of gradient evaluations used by each algorithm. The ESS of 
a set of M correlated samples 9^'^^ with respect to some function f{9) is the number of 
independent draws from the target distribution p{6) that would give a Monte Carlo estimate 
of the mean under p of f{0) with the same level of precision as the estimate given by the mean 
of / for the correlated samples 9^'^ . That is, the ESS of a sample is a measure of how many 
independent samples a set of correlated samples is worth for the purposes of estimating the 
mean of some function; a more efficient sampler will give a larger ESS for less computation. 
We use the number of gradient evaluations performed by an algorithm as a proxy for the 
total amount of computation performed; in all of the models and distributions we tested the 
computational overhead of both HMC and NUTS is dominated by the cost of computing 
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gradients. Details of the method we use to estimate ESS are provided in appendix [A} In 
each experiment, we discarded the first 1000 samples as burn-in when estimating ESS. 
ESS is inherently a univariate statistic, but all of the distributions we test HMC and 



NUTS on are multivariate. Following Girolami and Calderhead (2011) we compute ESS 



separately for each dimension and report the minimum ESS across all dimensions, since we 
want our samplers to effectively explore all dimensions of the target distribution. For each 
dimension we compute ESS in terms of the variance of the estimator of that dimension's 
mean and second central moment (where the estimate of the mean used to compute the 
second central moment is taken from a separate long run of 50,000 iterations of NUTS with 
6 = 0.5), reporting whichever statistic has a lower effective sample size. We include the 
second central moment as well as the mean because for simulation lengths eL that hit a 
resonance of the target distribution HMC can produce samples that are anti-correlated. 
These samples yield low-variance estimators of parameter means, but very high-variance 
estimators of parameter variances, so computing ESS only in terms of the mean of 9 can be 
misleading. 

4.1 Models and Datasets 

To evaluate NUTS and HMC, we used the two algorithms to sample from four target distri- 
butions, one of which was synthetic and the other three of which are posterior distributions 
arising from real datasets. 

250-dimensional multivariate normal (MVN): In these experiments the target dis- 
tribution was a zero-mean 250-dimensional multivariate normal with known precision matrix 
A, i.e., 

p{6) ocexpi-^e^Ae}. (20) 

The matrix A was generated from a Wishart distribution with identity scale matrix and 
250 degrees of freedom. This yields a target distribution with many strong correlations. 
The same matrix A was used in all experiments. 

Bayesian logistic regression (LR): In these experiments the target distribution is the 
posterior of a Bayesian logistic regression model fit to the German credit dataset (available 



from the UCI repository (Frank and Asuncion 2010)). The target distribution is 



p{a, f3\x, y) oc p{y\x, a, P)p{a)p{/3) ^^^^ 
oc exp{- log(l + exp{-yi{a + Xi ■ (3}) - - ^(3 ■ (3}, 

where 24-dimensional vector of numerical predictors associated with a customer 

i, yi is —1 if customer i should be denied credit and 1 if that customer should receive 
credit, a is an intercept term, and /3 is a vector of 24 regression coefficients. All predictors 
are normalized to have zero mean and unit variance, a and each element of (3 are given 
weak zero-mean normal priors with variance o"^ = 100. The dataset contains predictor and 
response data for 1000 customers. 

Hierarchical Bayesian logistic regression (HLR): In these experiments the target 
distribution is again the posterior of a Bayesian logistic regression model fit to the German 
credit dataset, but this time the variance parameter in the prior on a and /3 is given an 
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Figure 3: Discrepancies between the realized average acceptance probability statistic h and 
its target 5 for the multivariate normal, logistic regression, hierarchical logistic 
regression, and stochastic volatility models. Each point's distance from the x- 
axis shows how effectively the dual averaging algorithm tuned the step size e for 
a single experiment. Leftmost plots show experiments run with NUTS, other 
plots show experiments run with HMC with a different setting of eL. 

exponential prior and estimated as well. Also, we expand the predictor vectors by including 
two-way interactions, resulting in (^2^) + 24 = 300-dimensional vectors of predictors x and 
a 300-dimensional vector of coefficients p. These elaborations on the model make for a 
more challenging problem; the posterior is in higher dimensions, and the variance term o"^ 
interacts strongly with the remaining 301 variables. The target distribution for this problem 
is 

p{a, 13, a'^\x, y) oc p{y\x, a, f3)p{l5\a'^)p{a\a'^)p{a'^) 

oc exp{- Y,. log(l + exp{-yjXi • - - • /3 - f logcr^ - Aa^}, 

(22) 

where N = 1000 is the number of customers and A is the rate parameter to the prior on 
(j^. We set A = 0.01, yielding a weak exponential prior distribution on whose mean and 
standard deviation are 100. 

Stochastic volatility (SV): In the final set of experiments the target distribution is the 
posterior of a relatively simple stochastic volatility model fit to 3000 days of returns from 
the S&P 500 index. The model assumes that the observed values of the index are generated 
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Figure 4: Plots of the convergence of e as a function of the number of iterations of NUTS 
with dual averaging with 5 = 0.65 applied to the multivariate normal (MVN), 
logistic regression (LR), hierarchical logistic regression (HLR), and stochastic 
volatility (SV) models. Each trace is from an independent run. The y-axis shows 
the value of e, divided by one of the final values of e so that the scale of the traces 
for each problem can be readily compared. 



by the following generative process: 

r ~ Exponential(lOO); v ~ Exponential(lOO); si ~ Exponential (100); 

logSi>i ~ Normal(logSi_i,r-i); ^"^^'"'"^^-^ ~ t„ (23) 

where Si>i refers to a scale parameter Si where i > 1. We integrate out the precision 
parameter r to speed mixing, leading to the 3001-dimensional target distribution 

p{s,u\y) « e-0-0i-e-0-0i^i(n?£f t,(s-i(logy, -logy,_i)))x 

(0.01 + 0.5 E?=2°(log - log s,_i)2)-^ . (24) 

4.2 Convergence of Dual Averaging 

Figure [3] plots the realized versus target values of the statistics h^^'~' and h^^"^^. The h 
statistics were computed from the 1000 post-burn-in samples. The dual averaging algorithm 
of section |3.2| usually does a good job of coercing the statistic h to its desired value 5. It 
performs somewhat worse for the stochastic volatility model, which we attribute to the 
longer burn-in period needed for this model; since it takes more samples to reach the 
stationary regime for the stochastic volatility model, the adaptation algorithm has less time 
to tune e to be appropriate for the stationary distribution. This is particularly true for 
HMC with small values of 6, since the overly high rejection rates caused by setting 6 too 
small lead to slower convergence. 

Figure |4] plots the convergence of the averaged iterates Cm as a function of the number of 
dual averaging updates for NUTS with 6 = 0.65. Except for the stochastic volatility model, 
which requires longer to burn in, e roughly converges within a few hundred iterations. 
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Figure 5: Histograms of the trajectory lengths generated by NUTS with various accep- 
tance rate targets 5 for the multivariate normal (MVN), logistic regression (LR), 
hierarchical logistic regression (HLR), and stochastic volatility (SV) models. 



4.3 NUTS Trajectory Lengths 

Figure [5] shows histograms of the trajectory lengths generated by NUTS. Most of the tra- 
jectory lengths are integer powers of two, indicating that the U-turn criterion in equation 
[9] is usually satisfied only after a doubling is complete and not by one of the intermedi- 
ate subtrees generated during the doubling process. This behavior is desirable insofar as it 
means that we only occasionally have to throw out entire half-trajectories to satisfy detailed 
balance. 
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Figure 6: Effective sample size (ESS) as a function of b and (for HMC) simulation length 
eh for the multivariate normal, logistic regression, hierarchical logistic regression, 
and stochastic volatility models. Each point shows the ESS divided by the number 
of gradient evaluations for a separate experiment; lines denote the average of the 
points' y-values for a particular b. Leftmost plots are NUTS's performance, each 
other plot shows HMC's performance for a different setting of eL. 



The trajectory length (measured in number of states visited) grows as the acceptance 
rate target b grows, which is to be expected since a higher b will lead to a smaller step 
size e, which in turn will mean that more leapfrog steps are necessary before the trajectory 
doubles back on itself and satisfies equation [9] 
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Figure 7: Samples generated by random-walk Metropolis, Gibbs sampling, and NUTS. The plots 
compare 1,000 independent draws from a highly correlated 250-dimensional distribu- 
tion (right) with 1,000,000 samples (thinned to 1,000 samples for display) generated by 
random-walk Metropolis (left), 1,000,000 samples (thinned to 1,000 samples for display) 
generated by Gibbs sampling (second from left), and 1,000 samples generated by NUTS 
(second from right). Only the first two dimensions are shown here. 



4.4 Comparing the Efficiency of HMC and NUTS 

Figure [6] compares the efficiency of HMC (with various simulation lengths A ~ eL) and 
NUTS (which chooses simulation lengths automatically). The x-axis in each plot is the 
target 5 used by the dual averaging algorithm from section [3^ to automatically tune the step 
size e. The y-axis is the effective sample size (ESS) generated by each sampler, normalized by 
the number of gradient evaluations used in generating the samples. HMC's best performance 
seems to occur around 6 = 0.65, suggesting that this is indeed a reasonable default value 
for a variety of problems. NUTS's best performance seems to occur around 6 = 0.6, but 
does not seem to depend strongly on 6 within the range 5 G [0.45,0.65]. 6 = 0.6 therefore 
seems like a reasonable default value for NUTS. 

On the two logistic regression problems NUTS is able to produce effectively indepen- 
dent samples about as efficiently as HMC can. On the multivariate normal and stochastic 
volatility problems, NUTS with 6 = 0.6 outperforms HMC's best ESS by about a factor of 
three. 

As expected, HMC's performance degrades if an inappropriate simulation length is cho- 
sen. Across the four target distributions we tested, the best simulation lengths A for HMC 
varied by about a factor of 100, with the longest optimal A being 17.62 (for the multivari- 
ate normal) and the shortest optimal A being 0.17 (for the simple logistic regression). In 
practice, finding a good simulation length for HMC will usually require some number of 
preliminary runs. The results in Figure [6] suggest that NUTS can generate samples at least 
as efficiently as HMC, even discounting the cost of any preliminary runs needed to tune 
HMC's simulation length. 
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4.5 Qualitative Comparison of NUTS, Random- Walk Metropolis, and Gibbs 

In section [4^ we compared the efficiency of NUTS and HMC. In this section, we informally 
demonstrate the advantages of NUTS over the popular random-walk Metropolis (RWM) 
and Gibbs sampling algorithms. We ran NUTS, RWM, and Gibbs sampling on the 250- 
dimensional multivariate normal distribution described in section 14. 1[ NUTS was run with 
6 = 0.5 for 2,000 iterations, with the first 1,000 iterations being used as burn-in and to adapt 
e. This required about 1,000,000 gradient and likelihood evaluations in total. We ran RWM 
for 1,000,000 iterations with an isotropic normal proposal distribution whose variance was 



selected beforehand to produce the theoretically optimal acceptance rate of 0.234 (Gelman 



et al. , 1996). The cost per iteration of RWM is effectively identical to the cost per gradient 
evaluation of NUTS, and the two algorithms ran for about the same amount of time. We 
ran Gibbs sampling for 1,000,000 sweeps over the 250 parameters. This took longer to run 
than NUTS and RWM, since for the multivariate normal each Gibbs sweep costs more than 
a single gradient evaluation; we chose to nonetheless run the same number of Gibbs sweeps 
as RWM iterations, since for some other models Gibbs sweeps can be done more efficiently. 

Figure [7| visually compares independent samples (projected onto the first two dimen- 
sions) from the target distribution with samples generated by the three MCMC algorithms. 
RWM has barely begun to explore the space. Gibbs does better, but still has left parts 
of the space unexplored. NUTS, on the other hand, is able to generate many effectively 
independent samples. 

We use this simple example to visualize the relative performance of NUTS, Gibbs, 
and RWM on a moderately high-dimensional distribution exhibiting strong correlations. 
For the multivariate normal, Gibbs or RWM would of course work much better after an 
appropriate rotation of the parameter space. But finding and applying an appropriate 
rotation can be expensive when the number of parameters D gets large, and RWM and 
Gibbs both require 0{D^) operations per eff'ectively independent sample even under the 
highly optimistic assumption that a transformation can be found that renders all parameters 
i.i.d. and can be applied cheaply (e.g. in 0{D) rather than the usual 0{D^) cost of matrix- 
vector multiplication and the 0{D^) cost of matrix inversion). This is shown for RWM by 



Creutz (1988), and for Gibbs is the result of needing to apply a transformation requiring 
0{D) operations D times per Gibbs sweep. For complicated models, even more expensive 
transformations often cannot render the parameters sufficiently independent to make RWM 
and Gibbs run efficiently. NUTS, on the other hand, is able to efficiently sample from 
high- dimensional target distributions without needing to be tuned to the shape of those 
distributions. 



5. Discussion 

We have presented the No- U- Turn Sampler (NUTS), a variant of the powerful Hamilto- 
nian Monte Carlo (HMC) Markov chain Monte Carlo (MCMC) algorithm that eliminates 
HMC's dependence on a number-of-steps parameter L but retains (and in some cases im- 
proves upon) HMC's ability to generate eff'ectively independent samples efficiently. We also 
developed a method for automatically adapting the step size parameter e shared by NUTS 
and HMC via an adaptation of the dual averaging algorithm of Nesterov (2009), making 
it possible to run NUTS with no hand tuning at all. The dual averaging approach we 
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developed in this paper could also be applied to other MCMC algorithms in place of more 
traditional adaptive MCMC approaches based on the Robbins-Monro stochastic approxi- 



mation algorithm (Andrieu and Thoms, 2008; Robbins and Monro, 1951). 



In this paper we have only compared NUTS with the basic HMC algorithm, and not its 
extensions, several of which are reviewed by Neal (2011 ). We only considered simple kinetic 
energy functions of the form • r, but both NUTS and HMC can benefit from introducing 
a "mass" matrix M and using the kinetic energy function ^r'^ M~^r. If approximates 
the covariance matrix of then this kinetic energy function will reduce the negative 
impacts strong correlations and bad scaling have on the efficiency of both NUTS and HMC. 
Another extension of HMC introduced by Neal ( 1994 ) considers windows of proposed states 
rather than simply the state at the end of the trajectory to allow for larger step sizes without 
sacrificing acceptance rates (at the expense of introducing a window size parameter that 
must be tuned). The effectiveness of the windowed HMC algorithm suggests that NUTS's 
lack of a single accept/reject step may be responsible for some of its performance gains over 
vanilla HMC. 



Girolami and Calderhead (2011) recently introduced Riemannian Manifold Hamilto- 



nian Monte Carlo (RMHMC), a variant on HMC that simulates Hamiltonian dynamics in 
Riemannian rather than Euclidean spaces, effectively allowing for position-dependent mass 
matrices. Although the worst-case 0{D^) matrix inversion costs associated with this al- 
gorithm often make it expensive to apply in high dimensions, when these costs are not 
too onerous RMHMC's ability to adapt its kinetic energy function makes it very efficient. 
There are no technical obstacles that stand in the way of combining NUTS's ability to adapt 
its trajectory lengths with RMHMC's ability to adapt its mass matrices; exploring such a 
hybrid algorithm seems like a natural direction for future research. 

Like HMC, NUTS can only be used to resample unconstrained continuous- valued vari- 
ables with respect to which the target distribution is differentiable almost everywhere. HMC 
and NUTS can deal with simple constraints such as nonnegativity or restriction to the sim- 
plex by an appropriate change of variable, but discrete variables must either be summed out 
or handled by other algorithms such as Gibbs sampling. In models with discrete variables, 
NUTS's ability to automatically choose a trajectory length may make it more effective than 
HMC when discrete variables are present, since it is not tied to a single simulation length 
that may be appropriate for one setting of the discrete variables but not for others. 

Some models include hard constraints that are too complex to eliminate by a simple 
change of variables. Such models will have regions of the parameter space with posterior 
probability. When HMC encounters such a region, the best it can do is stop short and restart 



with a new momentum vector, wasting any work done before violating the constraints (Neal 



2011 ). By contrast, when NUTS encounters a 0-probability region it stops short and samples 



from the set of points visited up to that point, making at least some progress. 

NUTS with dual averaging makes it possible for Bayesian data analysts to obtain the 
efficiency of HMC without spending time and effort hand-tuning HMC's parameters. This 
is desirable even for those practitioners who have experience using and tuning HMC, but it 
is especially valuable for those who lack this experience. In particular, NUTS's ability to 
operate efficiently without user intervention makes it well suited for use in generic inference 
engines in the mold of BUGS (Gilks and Spiegelhalter , 1992), which until now have largely 
relied on much less efficient algorithms such as Gibbs sampling. We are currently devel- 
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oping an automatic Bayesian inference system called Stan, which uses NUTS as its core 
inference algorithm for continuous- valued parameters. Stan promises to be able to generate 
effectively independent samples from complex models' posteriors orders of magnitude faster 
than previous systems such as BUGS and JAGS. 

In summary, NUTS makes it possible to efficiently perform Bayesian posterior inference 
on a large class of complex, high-dimensional models with minimal human intervention. It is 
our hope that NUTS will allow researchers and data analysts to spend more time developing 
and testing models and less time worrying about how to fit those models to data. 
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Appendix A. Estimating Effective Sample Size 

For a function f{9), a target distribution p{9), and a Markov chain Monte Carlo (MCMC) 
sampler that produces a set of M correlated samples drawn from some distribution q{9^'^) 
such that g(6'™) = p{0^) for any m € {1, . . . , M}, the effective sample size (ESS) of 6^'-^ is 
the number of independent samples that would be needed to obtain a Monte Carlo estimate 
of the mean of / with equal variance to the MCMC estimate of the mean of /: 

Ess,j{e )-M i + 2T''-Hi — m7' 

f _ IE,[(/(gO - Ep[/(g)])(/(g*-^) - E,[/(g)])] 

where p( denotes the autocorrelation under (7 of / at lag s and Vp[x] denotes the variance 
of a random variable x under the distribution p(x). 

To estimate ESS, we first compute the following estimate of the autocorrelation spectrum 
for the function f{0): 



M 

' m=s+l 

where the estimates fi / and a'j of the mean and variance of the function / are computed with 
high precision from a separated 50,000-sample run of NUTS with 6 = 0.5. We do not take 
these estimates from the chain whose autocorrelations we are trying to estimate — doing 
so can lead to serious underestimates of the level of autocorrelation (and thus a serious 
overestimate of the number of effective samples) if the chain has not yet converged or has 
not yet generated a fair number of effectively independent samples. 

Any estimator of pi is necessarily noisy for large lags s, so using the naive estimator 
ESSo f(0^'^^) = ^.r ^, ——f will yield bad results. Instead, we truncate the sum over 
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the autocorrelations when the autocorrelations first dip below 0.05, yielding the estimator 
ESS,j{e'-'^') = -; Mf*°ff = min. s.t. p{ < 0.05. (27) 

We found that this method for estimating ESS gave more reliable confidence intervals 



for MCMC estimators than the autoregressive approach used by CODA (Plummer et al 



2006). (The more accurate estimator comes at the expense of needing to compute a costly 
high-quality estimate of the true mean and variance of the target distribution.) The 0.05 
cutoff is somewhat arbitrary; in our experiments we did not find the results to be very 
sensitive to the precise value of this cutoff. 
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